opts_knit$set(progress = TRUE, verbose = TRUE, root.dir = "~/Documents/git/ash/paper/Rcode")
require(ashr)

## Loading required package: ashr
## Loading required package: Rcpp

require(qvalue)

## Loading required package: qvalue

require(fdrtool)

## Loading required package: fdrtool

require(mixfdr)

## Loading required package: mixfdr

require(locfdr)

## Loading required package: locfdr
## Loading required package: splines

require(ggplot2)

## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## 
## The following object is masked from 'package:qvalue':
## 
##     qplot

load("sim1_out.RData")

#look into methods for estimating pi0 conservatively

# par(mfcol=c(3,3))
plot_ecdf = function(sims) {
    for (i in 1:length(sims$beta)) {
        plot(ecdf(sims$beta[[i]]), xlim = c(-6, 6), main = paste0("iteration ", 
            i))
        x = seq(-6, 6, length = 1000)
        lines(cdf.ash(sims$betahat.ash.n[[i]], x), col = 2, lwd = 2)
        lines(cdf.ash(sims$betahat.ash.u[[i]], x), col = 3, lwd = 2)
        lines(cdf.ash(sims$betahat.ash.true[[i]], x), col = 4, lwd = 2)
    }
}
plot_ecdf(simres1)

plot_ecdf(simres2)

Figure to show that estimated betahats are not so different

plot(betahat.ash.u[[1]]$PosteriorMean, betahat.ash.n[[1]]$PosteriorMean)

## Error: object 'betahat.ash.u' not found

abline(a = 0, b = 1, lwd = 2, col = 2)

## Error: plot.new has not been called yet

QUestion: is the null-biased prior maybe a little too conservative? Answer: log likelihoods don't suggest they are

# hh.ashtrue = hh.ashz hh.ashtrue$fitted.g$pi =
# c(2/3,1/15,1/15,1/15,1/15,1/15) hh.ashtrue$fitted.g$mean = c(0,0,0,0,0,0)
# hh.ashtrue$fitted.g$sd = sqrt(c(0,1,0.2,0.4,0.8,3))

# loglik(hh.ashtrue,betahat,sebetahat) loglik(hh.ashz,betahat,sebetahat)